      subroutine parcelc
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M, ONLY:
      use cophys_com_M
      use ctemp_com_M, ONLY: fmf, fef, fvf, fbxf, fbyf, fbzf, twx
      use blcom_com_M, ONLY: iphd2, wate, qdnc, iphead, ijkcell, ijkvtx, &
         ijkctmp, pxi, peta, pzta, link, qpar, periodicx
      use boundary
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
!
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer , dimension(8) :: iv
      integer , dimension(4) :: lioinput, lioprint
      integer , dimension(27) :: ijkstep
      integer :: l, n, ntmp, ntmp_nxt, np, ll, index, ijkl, ijkr, ntot
      real(double), dimension(1) :: ixi1, ixi2, ieta1, ieta2
      real(double) :: ixi4, ieta4, ixi5, ieta5
      real(double), dimension(8) :: wght
      real(double), dimension(100) :: fpxf, fpyf, fpzf, fpxft, fpyft, fpzft, &
         ucm, vcm, wcm, cm
      real(double), dimension(3,20) :: wsin, wcos
      real(double), dimension(3,12,20) :: tsi, rot
      real(double) :: nu, qptotl, the, zeta, wi, wim, wip, wj, wjm, wjp, wk, &
         wkm, wkp, qctotl, qtot2
      logical :: expnd, nomore
      character :: name*80
!-----------------------------------------------
!
!
!
!      a routine to interpolate particle data to the grid
!
      call celstep (iwid, jwid, kwid, ijkstep)
!
      do k = 1, kbp2
         do j = 1, jbp2
!dir$ ivdep
            if (ibp2 > 0) then
!
               qdnc(1+(j-1)*jwid+(k-1)*kwid:(ibp2-1)*iwid+1+(j-1)*jwid+(k-1)*&
                  kwid:iwid) = 0.0
               iphd2(1+(j-1)*jwid+(k-1)*kwid:(ibp2-1)*iwid+1+(j-1)*jwid+(k-1)*&
                  kwid:iwid) = 0
!
               ijk = (ibp2 - 1)*iwid + 1 + (j - 1)*jwid + (k - 1)*kwid
            endif
         end do
      end do
!
!dir$ ivdep
!
      wate(ijkvtx(:nvtx),:27) = 0.0
!
      ijkctmp(:ncells) = ijkcell(:ncells)
!
      ntmp = ncells
!
      qptotl = 0.0
!
    1 continue
      nomore = .TRUE.
!
!dir$ ivdep
!
      ntmp_nxt = 0
      do n = 1, ntmp
         if (iphead(ijkctmp(n)) == 0) cycle
         ntmp_nxt = ntmp_nxt + 1
         ijkctmp(ntmp_nxt) = ijkctmp(n)
         nomore = .FALSE.
      end do
!
      ntmp = ntmp_nxt
!
      if (.not.nomore) then
!
!dir$ ivdep
         do n = 1, ntmp
!
            ijk = ijkctmp(n)
            np = iphead(ijk)
!
!
!
            k = int(pzta(np))
            j = int(peta(np))
            i = int(pxi(np))
!
            the = pxi(np) - i - 0.5
            zeta = peta(np) - j - 0.5
            nu = pzta(np) - k - 0.5
!
            wi = 0.75 - the**2
            wim = 0.5*(0.5 - the)**2
            wip = 0.5*(0.5 + the)**2
!
            wj = 0.75 - zeta**2
            wjm = 0.5*(0.5 - zeta)**2
            wjp = 0.5*(0.5 + zeta)**2
!
            wk = 0.75 - nu**2
            wkm = 0.5*(0.5 - nu)**2
            wkp = 0.5*(0.5 + nu)**2
!
!     k-plane
!
            wate(ijk,1) = wi*wj*wk
            wate(ijk,2) = wip*wj*wk
            wate(ijk,3) = wip*wjp*wk
            wate(ijk,4) = wi*wjp*wk
            wate(ijk,5) = wim*wjp*wk
            wate(ijk,6) = wim*wj*wk
            wate(ijk,7) = wim*wjm*wk
            wate(ijk,8) = wi*wjm*wk
            wate(ijk,9) = wip*wjm*wk
!
!     k-1 - plane
!
            wate(ijk,10) = wi*wj*wkm
            wate(ijk,11) = wip*wj*wkm
            wate(ijk,12) = wip*wjp*wkm
            wate(ijk,13) = wi*wjp*wkm
            wate(ijk,14) = wim*wjp*wkm
            wate(ijk,15) = wim*wj*wkm
            wate(ijk,16) = wim*wjm*wkm
            wate(ijk,17) = wi*wjm*wkm
            wate(ijk,18) = wip*wjm*wkm
!
!     k+1 - plane
!
            wate(ijk,19) = wi*wj*wkp
            wate(ijk,20) = wip*wj*wkp
            wate(ijk,21) = wip*wjp*wkp
            wate(ijk,22) = wi*wjp*wkp
            wate(ijk,23) = wim*wjp*wkp
            wate(ijk,24) = wim*wj*wkp
            wate(ijk,25) = wim*wjm*wkp
            wate(ijk,26) = wi*wjm*wkp
            wate(ijk,27) = wip*wjm*wkp
!
         end do

!
!
!     calculate the charge density
!
!dir$ ivdep
!
            do n = 1, ntmp
!
               ijk = ijkctmp(n)
!
               qdnc(ijk + ijkstep(1:27)) = qdnc(ijk + ijkstep(1:27)) +&
                          qpar(iphead(ijk))*wate(ijk,1:27)
!
            end do
!
!
!
!
!     **************** finish boundary conditions **********************
!
         do n = 1, ntmp
            ijk = ijkctmp(n)
            np = iphead(ijk)
            if (np <= 0) cycle
            iphead(ijk) = link(np)
            link(np) = iphd2(ijk)
            iphd2(ijk) = np
!
         end do
      endif
!
      if (.not.nomore) go to 1
!
      do n = 1, ncells
!
         ijk = ijkcell(n)
!
         iphead(ijk) = iphd2(ijk)
!
         iphd2(ijk) = 0
!
      end do
!
	  if(periodicy) then
!
!     periodic in j
!
      do i = 1, ibp2
!dir$ ivdep
         do k = 1, kbp2
!
            ijkl = 1 + (i - 1)*iwid + (k - 1)*kwid
            ijkr = 1 + (i - 1)*iwid + jbp1*jwid + (k - 1)*kwid
!
            qdnc(ijkr-jwid) = qdnc(ijkr-jwid) + qdnc(ijkl)
            qdnc(ijkl+jwid) = qdnc(ijkl+jwid) + qdnc(ijkr)
!
         end do
      end do
!
	  else
!
!	  reflect in j
!
      do i = 1, ibp2
!dir$ ivdep
         do k = 1, kbp2
!
            ijkl = 1 + (i - 1)*iwid + (k - 1)*kwid
            ijkr = 1 + (i - 1)*iwid + jbp1*jwid + (k - 1)*kwid
!
            qdnc(ijkr-jwid) = qdnc(ijkr-jwid) + qdnc(ijkr)
            qdnc(ijkl+jwid) = qdnc(ijkl+jwid) + qdnc(ijkl)
!
         end do
      end do
!
      end if
!
      if (periodicx) then
!
!     periodic in i
!
!
         do j = 1, jbp2
!dir$ ivdep
            do k = 1, kbp2
!
               ijkl = 1 + (j - 1)*jwid + (k - 1)*kwid
               ijkr = 1 + ibp1*iwid + (j - 1)*jwid + (k - 1)*kwid
!
               qdnc(ijkr-iwid) = qdnc(ijkr-iwid) + qdnc(ijkl)
               qdnc(ijkl+iwid) = qdnc(ijkl+iwid) + qdnc(ijkr)
!
            end do
         end do
!
      else
!
!     reflect in i
!
!
         do j = 1, jbp2
!dir$ ivdep
            do k = 1, kbp2
!
               ijkl = 1 + (j - 1)*jwid + (k - 1)*kwid
               ijkr = 1 + ibp1*iwid + (j - 1)*jwid + (k - 1)*kwid
!
               qdnc(ijkr-iwid) = qdnc(ijkr-iwid) + qdnc(ijkr)
               qdnc(ijkl+iwid) = qdnc(ijkl+iwid) + qdnc(ijkl)
!
            end do
         end do
!
      endif

      if (periodicz) then
!
!     periodic in k
!
      do i = 1, ibp2
!dir$ ivdep
         do j = 1, jbp2
!
            ijkl = 1 + (i - 1)*iwid + (j - 1)*jwid
            ijkr = 1 + (i - 1)*iwid + (j - 1)*jwid + kbp1*kwid
!
            qdnc(ijkr-kwid) = qdnc(ijkr-kwid) + qdnc(ijkl)
            qdnc(ijkl+kwid) = qdnc(ijkl+kwid) + qdnc(ijkr)
!
         end do
      end do
      else
!
!     reflect in k
!
      do i = 1, ibp2
!dir$ ivdep
         do j = 1, jbp2
!
            ijkl = 1 + (i - 1)*iwid + (j - 1)*jwid
            ijkr = 1 + (i - 1)*iwid + (j - 1)*jwid + kbp1*kwid
!
            qdnc(ijkl+kwid) = qdnc(ijkl+kwid) + qdnc(ijkl)
            qdnc(ijkr-kwid) = qdnc(ijkr-kwid) + qdnc(ijkr)
!
         end do
      end do

      end if
!
!     diagnostics of some importance
!
      return
!
!      remove the `return` above to put the diagnostics back
!
      qctotl = 0.0
      qtot2 = 0.
      ntot = 0
!
      do k = 2, kbp1
         do j = 2, jbp1
            do i = 2, ibp1
!
               ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid
!
               qctotl = qctotl + qdnc(ijk)
!
               np = iphead(ijk)
  666          continue
               if (np <= 0) cycle
               qtot2 = qtot2 + qpar(np)
               ntot = ntot + 1
               np = link(np)
               go to 666
            end do
         end do
      end do
!
      write (*, *) 'parcelc: qctotl=', qctotl
      write (*, *) 'parcelc: qtot2=', qtot2
      write (*, *) 'parcelc: ntot=', ntot
      return
      end subroutine parcelc
